tf/QS/f C£~ /7?,a^3 


NASA Contractor Report 178033 

ICASE REPORT NO. 85-58 nasa-cr-178033 

1 9860009522 


ICASE 


FACTORIZATION AND REDUCTION METHODS 

FOR OPTIMAL CONTROL OF DISTRIBUTED PARAMETER SYSTEMS 


John A. Burns 
Robert K. Powers 


Contract Nos. NAS1-17070, NAS1-18107 
November 1985 


INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 


(VJASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 




FACTORIZATION AND REDUCTION METHODS 


FOR OPTIMAL CONTROL OF DISTRIBUTED PARAMETER SYSTEMS 


John A. Burns 

Virginia Polytechnic Institute and State University 

•k-k 

Robert K. Powers 

Institute for Computer Applications in Science and Engineering 


ABSTRACT 

A Chandrasekhar-type factorization method is applied to the linear- 
quadratic optimal control problem for distributed parameter systems. An 
aeroelastic control problem is used as a model example to demonstrate that if 
computationally efficient algorithms, such as those of Chandrasekhar-type, are 
combined with the special structure often available to a particular problem, 
then an abstract approximation theory developed for distributed parameter 
control theory becomes a viable method of solution. A numerical scheme based 
on averaging approximations is applied to hereditary control problems. 
Numerical examples are given. 
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I. INTRODUCTION 


The problem of developing computational algorithms for distributed 
parameter control systems has been the subject of a large number of recent 
papers. These articles run the gamut from very abstract papers dealing with 
general approximation theory for infinite dimensional systems to explicit 
numerical algorithms derived expressly for a particular type of system and/or 
application. In any investigation of approximation schemes for such systems, 
we believe that it is important to keep in mind the ultimate goal of the 
approximation. For example, an approximation scheme that produces excellent 
results if used for parameter estimation might be inappropriate for computing 
feedback gain operators. The point at which approximations are introduced 
into the analysis is something that also varies with the problem and with the 
particular approach used to analyze the algorithm. 

In this paper we concentrate on the linear quadratic optimal control 
problem for certain distributed parameter systems. We employ the 
approximation theory developed by Gibson [25] to formulate and analyze fast 
computational algorithms for approximating optimal feedback gain operators. 
These methods are based on factorization schemes of Chandrasekhar type. 

A primary objective of this paper is to show that if one combines a 
computationally efficient algorithm such as a Chandrasekhar type method with 
the special structure of often available in "real problems," then many 
practical problems can be attacked using distributed parameter control theory 
and sound computational techniques. As noted in Cast! and Ljung [15], this 
has been considered a major stumbling block between theoretical results for 
distributed parameter systems and their application to practical problems. 
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Our approach is based on an approximation theory specifically developed 
for approximating infinite dimensional control systems. Although the paper 
does contain some theoretical results on the existence, uniqueness, and 
smoothness of strong solutions to Chandrasekhar integral equations, we feel 
that the most significant aspects of this paper lie in the numerical results 
and applications. However, such theoretical results (especially smoothness) 
play an important role in the analysis of numerical schemes for direct 
integration of these equations (see Sorine [42]). These results have also 
been used to establish differentiability of strong solutions to operator 
Riccati differential equations (see [27], [44]). 

Previous authors have studied infinite dimensional Chandrasekhar equations 
in connection with quadratic control problems [8], [15], [16], [44], and many 
of these authors have suggested that such algorithms should lead to very 
efficient computational techniques. However, except for a few examples [12], 
[15], [16], [40], very few numerical results have appeared in the open 
literature that substantiate these claims. Therefore, we have included a 
number of numerical examples to illustrate the computational aspects of these 
algorithms. 

In Section 2 we present an aeroelastic control problem and briefly outline 
the derivation of the model. The control problem is then formulated as an 
infinite dimensional linear quadratic control problem. We introduce this 
model problem in order to provide an example to motivate our work and to check 
our computational algorithms. Section 3 is devoted to the development of the 
Chandrasekhar equations and computational algorithms for a general distributed 
parameter control problem. This section contains the statement of the major 
theoretical results. Proofs of these results are given in the Appendix. 
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In Section 4 we apply the general ideas developed in Section 3 to control 
problems governed by retarded functional differential equations. We 
restricted our attention to these problems for two reasons. First, except for 
some minor extensions, Gibson's approximation results [25] can be directly 
applied to obtain convergence of the Chandrasekhar algorithm. Secondly, we 
have considerable numerical experience with this class of problems, which 
allows us to compare the efficiency of a number of different numerical 
algorithms. We begin Section 4 by summarizing the basic results (see [2], 
[3], [25]) concerning the averaging approximation scheme for control problems 
governed by retarded functional differential equations. We also consider a 
"reduced averaging scheme" that takes advantage of the special structure that 
occurs in many problems and which often leads to considerable computational 
savings. Finally, we combine the Chandrasekhar algorithm with the reduced 
approximation scheme to produce a convergent and computationally efficient 
algorithm for approximating gain operators. The infinite time problem is also 
discussed. 

In Section 5 we present a number of numerical examples to illustrate the 
computational aspects of the theoretical results. These examples illustrate 
the potential applicability of these factorization methods for very large- 
scale control problems. As noted above proofs will be given in the Appendix. 

The notation used in this paper is fairly standard. Given two real 
Hilbert spaces X and Y, L(X,Y) shall denote the space of all bounded 
linear operators BtX + Y with the usual operator norm UBII. The 
innerproduct and the norm of the space X shall be denoted by <.,.> and 

A. 

H»ll x , respectively; subscripts will be dropped if it is clear from the 
context which space is intended. The adjoint of a closed linear operator 
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A:P(A)CX -»• Y shall be denoted by A*, The resolvent set of A is denoted 
by p(A). Finally, for an interval (a,b) of the real line, L 2 (a,b;X) 
denotes the set of L 2 -integrable functions whose values lie in the Hilbert 
space X. If the Hilbert space is understood from the context, L 2 (a,b) shall 
be used. 


II. AN AEROELASTIC CONTROL PROBLEM 

In order to provide the reader with some concrete example and to motivate 
our work, we present a brief description of a problem that we shall use to 
test the computational algorithms. A more detailed derivation of the model 
may be found in [10]. 

Consider the problem of controlling the aeroelastic structure (i.e., the 
typical section) shown in Figure 2.1. The airfoil is placed in a flow field 
with undisturbed stream velocity U and allowed to plunge and pitch in the 
f low. 



Figure 2.1 


Let h(t) denote the plunge and a(t) the pitch of the airfoil at time 


t. The equations of motion can be written in the form 


( 2 . 1 ) 


M y(t) + K y(t) = F(t) , 
s s 


T 

where y(t) = (h(t), a(t)) and F contains the aerodynamic and applied 
loads on the airfoil. In particular 


( 2 . 2 ) 


F(t) = 


L(t) + u(t) 




where L(t) and M (t) are the aerodynamic loads corresponding to the total 
wing lift per unit depth and total moment about the 1/4 chord per unit depth, 
respectively. For the airfoil considered here, it follows that (see [9], 
[10], [46]) 

(2.3) M a (t) = TT P b 3 [h(t)/2 + 3ba(t)/8 + Ua(t)] 
and 

(2.4) L(t) = npb^[h(t) + ba(t) + Ua(t)] + (2upUb)D(t) , 


where p is the density of air, b is the semichord length and D(t) is the 
"Duhamel integral." In particular, D(t) is given by 


(2.5) 



s))Q(s)ds 


where $(Ut/b) is the Wagner function (see [9]) and 
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In order to obtain a state space model that is suitable for control, it is 
necessary to provide a useful representation for D(t). Clearly, the key to 
this problem is the Wagner function $. It is possible to show that (see 
[ 11 ]) 

(2.7) D(t) = Q(t) - g(t) 

where g(t) is the output to a hereditary system (i.e., a functional 
differential equation) with input $(t). In particular, it follows that 

(2.8) *(Ut/b) = 1 - W(t) 
and 

(2.9) W(t) = Ce At B 

is the weighting pattern for a hereditary control system. Almost all 
approaches to modeling aeroelastic structures can be reduced to some scheme 
for approximating W(t). For example, R. T. Jones (see [30]) used the two- 
term exponential function 

-Bit -B 9 t 

(2.10) W(t) = e + a 2 e 

with ctj = 0.165, ct 2 ~ 0.335, = 0.0455 and ^ = 0.3. If W(t) is 

substituted into (2.5), then D(t) is approximated by 

(2.11) D(t) = Q(t) - g(t) = Q(t) - [ci 1 x^t) + a 2 x 2 (t)], 

where x^(t), x^(t) satisfy the second-order ordinary differential equation 
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with system parameters gj, q 2 , cj, and delay r > 0 (normally these 
parameters must be identified using a parameter estimation scheme). The 
problem of estimating the time delay r has been considered in other papers 
[4], [10]. The delay is clearly dependent on the chord length b and the 



undisturbed stream velocity U. We have found that for problems considered 
below a reasonable asymptotic "first estimate" of r is r = ^ where k = 
10. In this model, x d (t) represents a generalized aerodynamic "lag state." 
To complete the model we augment equations (2.1) - (2.4) with equation (2.14) 
and replace D(t) with the approximation D d (t) defined by (2.15). Let 
x(t) denote the five-dimensional vector 

(2.16) x(t) = (fi(t), a(t), h(t), a(t), x d (t)) T . 

The basic model becomes 

(2.17) Ex(t) = F q x(t) + F 1 x(t - r) + Gu(t) 
with initial data 

(2.18) x(0) = pe® 5 ; x d (s) = <f>(s)eL 2 (-r,0; ]R) 
and output 

(2.19) y ( t ) = Cx(t) € IBP. 

The control problem we consider is to find u :[0,T] K* that minimizes 

T 

(2.20) J(u) = J [ II y ( s ) II ^ + Ru^(s)]ds 

0 

where R > 0 and y(t) is the output to the delay-differential system (2.17) 
- (2.19). 



It is important to note that only the aerodynamic lag state x d (t) 
involves a time delay (i.e., the first four columns of F^ are zero). This 
observation will be important when we consider the computational algorithms in 
Section 5 below. The matrices E, Fq, Fp G, and C are given in the 
technical report [10]. Parameters needed to complete the modeling process can 
be obtained by applying the parameter identification methods presented in [10] 
to experimental wind tunnel data. A specific example of this process is given 
in [10]. 

There are a number of state space formulations for this problem. The 
approach we take here is a slight variation of the "standard" formulation 
given in references [9], [24], [46], Let Aq = E 1 Fq, A^ = E 1 F ^ , 

B = E * G and note that A^ has the form 


( 2 . 21 ) 




> 


where A d is a 5x1 (i.e., column) vector. We choose H - IB? x L2(-r,0;lR) 
as our state space and let (n»<t>( •)) = (ri^ , r^, 4>( * ) ) denote a 

typical element in H. Define the linear operator A on H 


(2.22) V(A) = {(n, <j>(O)|<|>(-)eW 1 ’ 2 (-r,0;K), <(,(0) = n 5 } 
by 

(2.23) A(n, 4.(0) = (A 0 n + A d 4>(-r), $(•)). 
Moreover, we define B:H + H and V:H IP by 


(2.24) 


B u = (Bu,9) 
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and 

(2.25) V(n, <(>(•)) = Cn, 

respectively. The delay-differential system (2.17) - (2.19) can be realized 
(see [10]) as the system 

(2.26) z(t) = Az(t) + 8u(t) 

(2*27) z(0) = z Q = (n»4>( *)) 

(2.28) y(t) = Vz(t) . 

Moreover, the optimal control problem is equivalent to finding u*:[0,T] -* IR*- 
that minimizes 

T 

(2*29) J(u) = / [ H y ( s ) II 2 + Ru 2 (s)]ds, 

0 

where y(t) is the output to (2.26) - (2.28). 

The problem defined by (2.26) - (2.29) will be used to test the numerical 
schemes described in Section 4. However, it also serves to motivate the 
theoretical developments presented below. 


III. THE CHANDRASEKHAR ALGORITHM 

In this section our attention is focused on the time-invariant infinite 
dimensional linear regulator problem. For a study of control problems 
governed by a general evolution equation we refer the interested reader to 
Curtain and Pritchard [18], Gibson [26], and Datko [19]. The evolution 
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processes usually arise from control of systems governed by partial 
differential equations (see Lions [35] and Lukes and Russell [36]), or 
functional differential equations (see Delfour and Mitter [23] and Manitius 
[38]). The approach we follow is similar to the development of the optimal 
control problem in Gibson [25]. The importance and potential usefulness of 
Chandrasekhar equations have been known for some time. There are a number of 
papers that discuss the application of these equations to finite dimensional 
linear quadratic optimal control problems (i.e., see [13], [14], [31], 
[22]). However, the derivation of these equations relied upon being able to 
twice differentiate the solution to a matrix Riccati differential equation. 
In an infinite dimensional setting there are a number of questions concerning 
the existence of these derivations (in a strong sense) that limit the 
usefulness of this approach for infinite dimensional systems. The derivation 
of the infinite dimensional equations in [45] by Tung is purely formal. The 
differentiations used in his derivation were not justified. In [8] and [15] 
the authors used a Lions-type setting and derived a set of Chandrasekhar 
differential equations satisfied in a distributional sense. We intend to give 
an alternate derivation which will lead to a set of Chandrasekhar integral 
equations that have unique strong solutions. Moreover, under fairly weak 
assumptions it can be shown that the solution to these integral equations is 
strongly differentiable (a result that can be applied to establish smoothness 
properties of solutions to Riccati operator differential equations [27]). 

In [44], Sorine established that the gain operator satisfies a set of 
Chandrasekhar equations if the underlying semigroup of the system is 
analytic. For systems governed by hyperbolic PDE's or differential-delay 
equations (such as the aeroelastic system above) the associated semigroups are 
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not analytic, and hence Sorine's results do not apply. In this section we 
give a set of Chandrasekhar equations in an integral form for a large class of 
optimal control problems governed by a general distributed parameter model. 
The results are established by an approximation technique and leads to 
computationally feasible methods. Because the proofs are technical and 
provide little insight into the structure of the problem, we have placed them 
in an appendix. 

Let H and U be real Hilbert spaces. Throughout this section 
T(t):H + H, t ^ 0, will denote a strongly continuous Cg-seraigroup of bounded 
linear operators with infinitesimal generator A. We shall always assume that 
the following basic hypothesis holds. 

The operators 8:U + H, Q :H * H and R:U -*• U are continuous 
linear operators. 

The operators Q and R are self-adjoint and non-negative and 
R satisfies IIRII ^ m > 0. 

The operator Q can be factored into the form Q = V* V where 
V:tf A is a bounded linear operator and A is a Hilbert space. 

The infinite dimensional linear quadratic (LQ) optimal control problem 
is to find a u €L 2 ( 0 ,t^;U) which minimizes 

C f 

(3»1) J(z,u) = f [<Qz(s) , z(s)> + <Ru(s) ,u(s)>1ds 

0 


H 0 ) i) 
ii) 
iii) 


where z(t) is defined by 
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t 

(3.2) z(t) = T(t)z + / BT(t - s) u(s)da 

0 

for 0 £ t £ and zeH* Although our discussion is mainly restricted to 
the control problem on finite intervals, we will have occasion to discuss its 
relationship to the infinite interval control problem. In this case the cost 
functional (3.1) becomes 

00 

(3.3) J (z,u) = / [<Qz(s),z(s)> + <Ru(s) ,u(s)>]ds. 

0 

Under the assumptions given above on Q and R, it is known that (see 

ic 

[18], [25] for example) there exists a unique u £L 2 (0,t f ;U), (respectively 
L 2 (0,t;U) for each t > 0) which minimizes J(z,u) (respectively 
J^Zqju)). The following characterization of this optimal control may be 
found in Gibson [25]. This characterization is discussed in detail in order 
to develop the notation we shall need for the derivation of the Chandrasekhar 
equations. 

For 0 < s < t, define the Hilbert spaces H and U by 

“ L 2 (s,t f ;H) and U g = L 2 (s,t f ;U), respectively. Let the operators 

T €L(H,H ), T €L(H ,H), and F €L(H ,H) be defined by 
o s s s s s s 

(T z)(t) = T(t - s)z, zeH, 

b 

t 

(Tg 4>)(t) = / T(t - r»)<l)(n)dri , 0 < s < t < t f , 4 »cH s , 

r s t = (T s <(>)( t f ) , 


and 
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respectively. Straightforward calculations imply that T is given by 

s 


CT t^)< t) = / T (p - t)i|)(n)dp, \peH 
t s 


and that F has the form 
s 


(F* z)(t) = T*(t, - t)z , z€H. 


These representations will be useful in the proofs of Lemmas 3.1 and 3.2 in 
the Appendix. 

For the optimal control problem on the interval [s,t f ] with 0 < s < t f , 
define the cost functional 

fc f 

J(s,z(s),u) = + / (<Qz(p),z(p)> + <Ru(n),u(n)>)dp, 

S 


where z(t) is defined by (3.2). In [25] Gibson has shown that the optimal 
control may be expressed in the form 

(3.4) u*(t) = -(r" 1 B*)(t)z*(s), 0 < s < t, 

A Ak 

for almost all t in [O-tf] where R and B are defined by 

r s s J 


(3.5) 

and 


R = R + B* T* QT 
s s s 


* * 

B s ■ 8 T s QT S > 


(3.6) 
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respectively. Note that R£l(U ,11 ) and B EL(H 9 U ). Since T(t)z - 

s s s s s 

T(t - s)z*(s) it follows directly from (3.2) and (3.4) that the optimal 
trajectory z (t) has the representation 


(3.7) z*(t) = T(t - s)z*(s) - / T(t - n)8(R" 1 8*)(p)z*(s)dp, 

S 5 

s 


for 0 _< s < t < t f . If S(t,s )eL(H,H) is defined by 


(3.8) 


S(t,s)z = T(t - s)z - / T(t - B*)(n)zd n , 

s 


then (3.7) becomes 


z (t) = S(t,s)z, 0 _< s <_ t _< t^. 


In [25], Gibson shows that S(t,s) is the bounded perturbation of T(t) 
— 1 

by -BR 8 n(t) where Il(t) is defined by 


(3.9) 


n(t)z = / T (p - t)QS(n,t)zdn, z etf, 
t 


Furthermore, he shows that the optimal control has the representation 


u*(t) = -R -1 B* n(t)z*(t) = -K(t)z*(t) . 


It has also been shown (see [18] and [25]) that Tl(t) is the unique self- 
adjoint solution to the two Riccati-type integral equations 


n(t)z = / T*(s - t) [Q - H(s)BR 1 8* H(s)]T(s - t)zds 

t 


(3.10) 
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and 
(3.11) 

A formal differentiation of either of these two equations yields a Riccati 
differential equation analogous to the case for the finite dimensional linear 
regulator problem. Such formal differentiations will be avoided in our 
derivation of the infinite dimensional Chandrasekhar equations. 

As in the finite dimensional case, under appropriate stabilizability and 
detectability requirements, the solution of the infinite time problem may be 
viewed as the limit of the solution to the finite time problem as t f ->■ «. 
Our only interest in this problem is to show that the Chandrasekhar equations 
may also lead to a computationally feasible method for computing solutions to 
these problems. Numerical examples will be given in Section 5 to illustrate 
this idea and to discuss some of its limitations. 

Since our main objective is the development of computational algorithms 
for numerically calculating the feedback gain operator for the optimal control 
problem (3.1) - (3.2), we shall approach the development of the Chandrasekhar 
equations via approximation theory. This theory is based on Gibson's work 
[25], [26], and our presentation makes use of a slight extension of his 
results. Again, we want to emphasize that these theoretical extensions are 
minor and that the real advantage of our approach lies in the computational 
savings that come from employing the problem structure and Chandrasekhar 
methods. However, this extension will enable us to show that the gain 
operator K(t) for the system (3.1) - (3.2) satisfies, in a strong sense, a 
set of Chandrasekhar integral equations. Moreover, examples will be given to 
show that approximation of the gain through the Chandrasekhar equations can be 


H(t)z 


S(s,t) [Q + n(s)BR _1 b" n(s) ]S(s,t)zds. 
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computationally more efficient than the usual method of solving a Riccati 
equation. 

Let {T N (t)} t>0 be a sequence of CQ-semigroups on H with corresponding 
infinitesimal generators {A N }. We also assume that B N eL(U,H) and Q N £L(H>H) 
define sequences of operators with each self-adjoint and non-negative. 

The following hypothesis is needed. 

Hj) i) There exist constants Tq, 1^, r 2 such that 

ii t n ( t ) ii £ r Q 

for all N and tc[0, t£ ] and 

Id I < Ij, II Q N II _< r 2 

for all N. 

ii) The semigroups T N (t) and [T N (t)j* converge strongly to T(t) 
and T (t), respectively and the convergence is uniform in t on 
[Ojt^J . 

iii) The operators B N , [B N ] and Q N converge strongly to 6, 8* 
and Q, respectively. 

It is important to note that condition H ^ - (ii) requires convergence of 
the semigroups and their adjoints. This is a critical assumption if one is 
interested in applying these results to non-self-ad joint problems. 

N A 

The Nth approximate LQ optimal control problem is to find [u ] (t) 


which minimizes 
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(3.12) J N ( Z N , u N ) = | [<Q^ z^(s),z^(s)> + <Ru N (s) ,u^(s)>]ds 

0 


where z^(t) is defined by 


(3.13) z N (t) = T N (t)z(0) + / T N (t - u N (s)ds 

0 

for 0 t < t . Under the above assumptions, Gibson [25] has shown that the 

- S 

N A 

optimal control [u ] converges in L 2 (0,t f ;U) to the optimal control u* 
for the original problem. In particular, he proved that the Riccati operators 
n N (t) converge strongly to n(t), and uniformly In t on compact t-intervals 
This yields strong convergence of the gains K N (t) = R _1 B N n N (t) to 
K(t) = R 8 n(t) . If the control space is finite dimensional, the convergence 
of K (t) is uniform in norm. 

In order to establish analogous results for Chandrasekhar equations, one 
needs the following technical lemmas, the proofs of which appear in the 
Appendix. 


LEMMA. 3.1: _If_ Hq and are satisfied, then l|B^[R^] ^(t)ll is 

uniformly bounded in N, s, and t for 0 s <_ t t^. 

Gibson [25] established that conditions Hq and H ^ are sufficient to 
ensure the strong convergence of S N (t,s) to S(t,s). 
slightly more general result. 


We make use of a 
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LEMMA 3.2: JL£ Hq and H ± hold, then for each z etf , 0 s ^ t ^ t^ 

N 

S (t,s)z -»■ S(t,s)z 

and 

[S N (t,s)] z + S*(t,s)z, 

where the convergence Is uniform In s and t. In addition , S^(t,s) and 
[S^(t,s)] are uniformly bounded in N, s, t for 0 _< s < t <_ t^. 

In order to derive an infinite dimensional version of the Chandrasekhar 
equations, it is necessary to make an additional assumption on the 
approximating sequences. This condition essentially requires that each of the 
approximating LQ control problems (3.12) - (3.13) have optimal gain operators 
that satisfy a form of the Chandrasekhar differential equations. 


H 2 ) i) There exists a sequence of approximating LQ optimal control 
problems (3.12) - (3.13) satisfying 
ii) The operators can be factored into 0^ - [V^] with 

+ V and [V^] •> V strongly. 

iii) The optimal gain operator K^(t) for the problem (3.12) - (3.13) 
is strongly continuously differentiable and satisfies 


(3.14) At)z - -R" 1 8 N *[L N (t)]* L N (t), K N (t f ) « 0 

(3.15) t N (t)z = -L N (t)[A N - B N K N (t)]z, L N (t f ) = V N 



- 20 - 


for zeH, where is the infintesimal generator of T^(t) and 

L N (t) - V N s"(t £; t). 

Condition H 2 ) - (iii) i s always satisfied if each A N is a bounded 

linear operator. In this case the semigroups T N (t) are differentiable, and 
the derivation of (3.14) - (3.15) proceeds as in the finite dimensional 

case. It should also be noted that conditions H j) - (iii) an d H 2 ) - (ii) 
are independent in that strong convergence of Q N to Q does not imply H 2 ) 
(■^i) and conversely that H 2 ) ~ (iii) is not sufficient to imply strong 
convergence of Q N to Q. 

We may now state the fundamental existence result. 


THEOREM 3.3: Suppose that conditions Hq, H ^ and hold. If K(t) 

denotes the optimal gain operator for problem (3.1) - (3.2), then K(t) 
satisfies the system of Chandrasekhar integral equations 


(3.16) 


K(t)z = / R 1 B* L*(s)L(s)zds 


(3.17) 


L(t)z = VT(t f - t)z - J L(s)8K(s)T(s - t)zds 

t 


f_or a ,ll z€H an ^ 0 S. t an ^ L(t)z = VS(t^,t)z. Moreover, the 

approximate gain operators K N (t) converge strongly to K(t). If U is 
finite dimensional, then 


lim UK At) - K(t)« = 0 


e convergence is uniform in t on [0,tj]. 


(3.19) 



- 21 - 


Theorem 3*3 does not imply that K(t) and L(t) are the only strongly 
continuous solutions to (3,16) - (3.17). The proof that these equations have 
unique strongly continuous solutions is non-trivial, but essential if we hope 
to use these equations as a basis for computational algorithms. 

THEOREM 3%4: Assume that Hq, H^ and H£ hold. Then K(t) and L(t) = 
VS(t£>t) are the unique strongly continuous solutions to the Chandrasekhar 
equations (3.16) - (3.17) . 

Some comments concerning the assumptions Hq, H^ and H 2 are in 
order. Although these conditions may at first glance seem severe, they are in 
fact the properties one would like to have in a scheme that is to be used to 
numerically approximate (3.1) - (3.2). In this case the convergence 
properties are assured and the Chandrasekhar integral equations may be 
approximated directly in order to obtain the gain operator. The averaging 
approximation scheme which is discussed in Section 4 satisfies these 
conditions when applied to linear regulator problems governed by delay- 
differential equations. For distributed parameter systems, Lukes and Russell 
[36] state conditions which also satisfy these hypotheses. Their 
approximations are essentially eigenfunction expansions and apply to problems 
whose infinitesimal generators have normal extensions (i.e., heat and wave 
equations with the appropriate boundary conditions). 

While our major concern is to illustrate that approximation of the gain 
operator through the Chandrasekhar equations may have computational advantages 
over approximation of the Riccati operator, it should be pointed out that the 
theory developed is not restricted to finite dimensional approximations. This 
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led Kazufumi Ito (see [27]) to note that if condition Hq was satisfied, then 
the Yosida approximations could be used to define approximating sequences that 
satisfy and H 2 • In particular, let A N = NA(NI - A) _1 (for N 
sufficiently large so that N€p(A>) and define B N =8, Q N = Q, and V N = V 
for each N. It can be shown that these approximating sequences satisfy H ^ 
and H 2 and hence one has the following theorem on the existence and 
uniqueness of Chandrasekhar integral equations. 

THEOREM 3.5: lf_ Hq is satisfied, then K(t) and L(t) = VS(t f ,t) are 
the unique strongly continuous solutions to the Chandrasekhar equations (3.16) 
~ (3.17). 


Although the Yosida approximates provide a set of approximating sequences 
that satisfy H 1 and H 2 , in general these approximations are not very useful 
for numerical schemes. For example, the A^ operators are infinite 
dimensional and not easily constructed. It is for this reason that we have 
stated the conditions that are needed to not only obtain existence and 
uniqueness but to also guarantee the convergence of the gain operators. 

We turn now to the application of these ideas to a particular class of 
distributed parameter systems. In particular, we shall concentrate on control 
systems governed by functional differential equations. 
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IV HEREDITARY CONTROL PROBLEMS 


In this section we consider the control of a linear hereditary system of 
the form 


(4.1) 


x(t) - Lx fc + Bq u(t), 0 t t^ 


x(o) = n, x Q = <j>, 


where n £ TEf 1 , ( — r , 0 ; I?) , u£L 2 ( 0 ,tj; iP) for each t^ < », B£L(lP, Ip) , 
and x = x(t + 9) for -r 0 0. The linear operator L has the form 

v 0 

Lx = l A x(t - h.) + / D(0)x(t + 0 )d0 , 

i=0 1 -r 

where each A^LClP.lP), 0 = h Q < <•••< h y = r, and D£L 2 (-r,0;L(TBf\lP)). 

Let Q be a real, symmetric, nonnegative n x n matrix, and R a real, 
symmetric, positive m x m matrix. The optimal control problem is to find 

•k 

u (t)£L 2 (0,t^; iP) which minimizes 

t f 

(4.2) J(u) / [<Qx(s) ,x(s)> — <Ru(s) ,u(s)>]ds 

0 

where x(t) is the solution to (4.1) corresponding to u(t). 

It is well-known that (see [3] for example) the system (4.1) has a unique 
solution for each u(t)£L 2 (0, iP) and (n,<j>)£lRx L 2 (-r,0;lP) = Z. In 
addition, the solution depends continuously on q, <j>, and u. Under the above 
stated conditions on Q and R there exists a unique u*€L 2 (0,t^; IP) 
which minimizes (4.2) subject to (4.1). 
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In recent years a standard technique for obtaining the solution to (4.1)- 

(4.2) has been to reformulate the system as an evolution equation on the 

Hilbert space Z. The solution to problem (4.1) - (4.2) is then approximated 
via the abstract formulation. Taking this approach, we will see that the 

dynamical system falls into the framework described in the previous sections 
and the optimal gain operator is thus characterized via the Chandrasekhar 

equations (3.16) - (3.17). 

For (n,<f>(*))£Z, define the operator T(t )£L(H,H) by 

(4.3) T(t)(n,<j>( •)) = (x(t),x t (»)) 

where x(t) is the solution to the homogeneous version of (4.1), i.e., 

u(t) = 0. It is well-known (see [3]) that T(t) is a Cg-semigroup on Z 

with infinitesimal generator A characterized by 

(4.4) 0(A) - {(n,*(O)|*(O€W 1 * 2 (-r,0;if l ), n = <f>(0) } 
and 

(4.5) A(ri, <}>(»)) = (L<(»,iJ(.)) 

for (n,<f>( •))££>( A). Moreover, the delay differential equation (4.1) is 
equivalent to the abstract evolution equation (AEE) 

t 

(4.6) z(t) = T(t)(ri,(f)( *)) + / T(t - s)Bu(s)ds 

0 

where Bu = (B^ u,0). In particular, the following result may be found in 


[ 4 ]. 
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THEOREM 4.1: Let (n,ij>)€Z be given. If x(t;u) is the solution of 

(4.1) for ueL 2 (0,t f ;lP), then z(t) defined by (4.6) satisfies z(t) = 

(x(t;u), x t (*;u)), t ^ 0. 


Define the operator Q€l(Z,Z) by 


Q = 


"Q 

0 


O' 

0 . 


where 0 denotes the zero operator on the appropriate spaces. The control 
problem (4.1) - (4.2) is equivalent to finding u (• )£L 2 (0, t^; iP) which 
minimizes 

tf 

(4.7) J = / [<Qz(s) ,z(s)> + <Ru(s) ,u(s)>]ds 

0 

where z(t) is defined by (4.6). Since the problem defined by (4.6) - (4.7) 
falls into the framework above, the optimal control has the representation 


(4.8) 


u*(t) = -R 1 8* n(t)z*(t) = — K(t)z*(t) 


where n(t) satisfies the Riccati integral equations (3.10) and (3.11) and 
z (t) is the optimal trajectory. Due to the special structure of the state 
space Z - IF? x L 2 (-r,0 jffP), n(t) may be expressed as the matrix 
of operators 


l~ A A 


n oo (t) 

n 01 

- n i0^ t ^ 

A 

n ll 


(t) 

(t) 


0 < t < t f 


(4.9) 


n(t) = 
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where there exists a real, nonnegative, symmetric n x n matrix 

A A 

n 00 (t)n = ^00^^* n io (t) a ^ oun< ied linear operator from 

and can be represented by 


n 00 ^ 

L 2 (-r,0;I?), 


[n 1Q (t)n](s) n io(t>s)n 


where n lQ (t,s) is an n x n matrix function, H Q1 (t) is the bounded linear 
operator from (-r,0; 2? 1 ) into ]FP given by 


(4.10) 


n 0 i(t)<{> = / n 10 (t;s)(j)(s)ds, 


-r 


10 ' 


^(-r.Ojlf 1 ), 


and n n (t) is a real, non-negative, self-adjoint operator on L2(-r,0; IF? 1 ) . 
Further properties of these operators may be found in [25], [33], and [35]. 

Combining (4.8) - (4.10), it follows that the optimal control may be written 
as 


(4.11) 


u*(t) = -R -1 Bj[n 00 (t)x*(t) + / IIj 0 (t,s)x*(s)ds] 

-r 


f° r 0 Therefore, the gain operator has the form 

0 

< 4 * 12 ) K(t)(n, $(•)) = K 0 (t)n + / K 1 (t,s)<j>(s)ds 

-r 

where Kg(t) = R RqqC t ) is an m x n matrix and 

— 1 T T 

K^Ct.s) = R B Q n 1Q (t,s) is an m x n matrix valued function. For each 
t >_ 0 the function s -*■ KjCtjS) is called the functional gain . 

It is important to point out that the Riccati operator maps an infinite 
dimensional space to itself while the gain operator maps an infinite 
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dimensional space to a finite dimensional space. Note that the gain is the 
actual operator needed to produce the optimal control. Thus it should be 
advantageous to approximate K(t) instead of n(t) since it involves 
approximations in one spatial variable only. This difference can be accounted 
for by examining equations (4.9) and (4.10). The infinite dimensionality in 

A A 

two spatial variables of H(t) is due to the operators n^(t) and H^o ( t ) • 

A 

However, the operator TI^(t) does not appear in the representation (4.11), 

A yv 

and the information in ll^t) is contained in the operator ^(t). The 
Chandrasekhar equations exploit this reduced dimensionality, and for the 
particular approximation scheme below it will be seen that this reduction can 
be quite significant. 

We turn now to a particular approximation scheme (the so-called AVE scheme 
[3], [4], [25]) for the simple delay differential equation 

x(t) = Aq x(t) + A^ x(t - r) + B q u(t) 

(4.13) 

x(0) = q, x(s) = ij)(s ) , -r _< s < 0. 

This scheme has been extensively studied by Banks and Burns [3] and Gibson 
[25], The treatment of a more general case which includes multiple discrete 
delays and a continuously delayed term may be found in [3]. 

For any positive integer N, partition [-r,0] into subintervals 
[tj, tj_^ ] for j = 1 ,2, • • • ,N, where t*? = - Let denote the 

characteristic function on [t^, tj_j) for j = 2,3, and the 

characteristic function on [-r/N,0]. Define the finite dimensional 
subspaces Z^c Z by 

z = { (b , <j>)ez j qe IP, * - ^ v. X *J, Vj e iP} 
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and the projections 


N 

P :Z Z 


N 


by 


An,*) - (*”, l ♦; x?) 

1=1 J 


where 


t N 

i-1 

*0 = n » " 7 / *(s)ds, for j = 1,2,' 

N 

j 

N N 

Let the operators A :Z + Z be defined by 


,N. 


.N, _ r . ,N N N n , N Ni N, 

A (u,*) - (A q * 0 + A l * N , - l (<(, - *.)x. )• 

j=l J J J 


If T N (t) denotes the semigroup generated by A N , then the approximation to 
the abstract formulation of (4.13) is 


z N (t) = T N (t)P N (n,<f>) + J T N (t - a)(B n u(a),0)da. 

0 U 


Since A^ 
problem 


is bounded for each N, this is equivalent to the initial value 


(4.14) 


z N (t) = A N z N (t) + (B q u(t),0) 
z N (0) = P N (q,*) 


Note that the operator A N is reduced by the Hilbert space Z N for 
each N. Since (Bgu(t) ,0)£Z^ for every N, (4.14) is a differential 
equation in a finite dimensional space. Upon choosing the appropriate basis, 
equation (4.14) has the representation (see [3]) 
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(4.15) 


•N,.* .N N, . . _N . . 
w (t) = A w (t) + B u(t) 

N, m w N N N> 

w (0) = col (4> q , 


where A N is the n(N + 1) x n(N +1) matrix 


(4.16) 


— I 
r 


N 

r 


0 “ I - “ I 

r r 


I is the n x n identity matrix, and is the n(N + 1) x m matrix 

B N = col(B 0 ,0,...,0). 

The Nth approximate control problem becomes: find the 

[u N ] €1^(0,^; 1^) which minimizes 

fc f 

(4.17) J N (w N (0),u N ) = / (<Q N w N (s),w N (s)> + <Ru N (s),u N (s)>)ds, 


where w^(t) is the solution to (4.15). The matrix is given by 


q n = 


Q o 

0 


• • • 


where 0 denotes the appropriate n x n zero matrix. 
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In [3] the authors proved that the AVE scheme satisfies a Trotter-Kato 
type approximation theorem. This established the strong convergence of the 
semigroups T N (t) to T(t) defined by (4.3). Recall that one of the 
requirements of H ^ is that the adjoint semigroups [T^(t)]* converge 

strongly to T (t). This result was obtained by Gibson [25] and is the key 
Inct used to prove Lemma 3.2. The special forms of and immediately 

satisfy and thus the solutions to the Chandrasekhar equations 

(3.14) - (3.15) of the AVE approximation converge to the solutions of the 
Chandrasekhar integral equations (3.16) - (3.17) in the strong sense. 

Before demonstrating that approximation of the gain operator via the 
Chandrasekhar equations has significant computational advantages over 

approximation of the gain operator via the Riccati equation, we introduce an 
additional reduction technique which takes advantage of a special structure 
that frequently occurs in hereditary systems. For these special systems, this 
technique will further reduce the number of differential equations that it is 
necessary to solve in order to approximate the gain operator, K(t). 

In applications of hereditary control problems, the delay does not always 
appear in each component of the state (see Examples 5.1, 5.2, and 5.3). In 
this type of system the initial data actually lie in a "smaller" state 
space. The AVE scheme as described projects each of the n components of the 
history portion of (x(t),x fc ) into an N-dimensional space. If only q < n 
of the components are delayed, it should be possible to project the history 
portion into a space of dimension q x N instead of a space of dimension 
n x N. The following discussion develops this idea. 

The reduction as discussed here was originally studied in Cliff and Burns 
[17] for spline approximations in the context of parameter identification, and 
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then for AVE in [10]. This reduction is a special case of the F-reduction 
discussed in [21] and [22]. 

For the remainder of this section it will be assumed that the matrix A^ 
defined in equation (4.13) has the form 



where Aj^ is p x q, A^ is q x q, and p + q = n. We shall write x(t) - 

col(w(t) ,y(t))€I? x T? and not distinguish between columns and rows for ease 

of notation. Also, we shall not distinguish between the space Z defined 

earlier in this section and IE? x 1$ x L^C-r,!);!?) x L^C-rjO; • 

Let Z^ denote the U reduced space 11 1? x x L£ (“r ,0 ; !&) . Define the 

projection operator P:Z + Z D by 

K 

P( H ^ j H2 > 4* y y ^2 ^ ~~ ^ D y * ^2 * ^2 ^ y 

and the injection I:Z + Z by 

K 

A 

^ (n 1 >^2 > ^2^ = \ > ^2 9 d y 9 $2 ^ * 

A 

where is the function in L^C-r ,0; I?) with constant value r^. For 

A and T(t) as defined by (4.3) - (4.5) define the reduced operator 

V«V CZ R* Z R b ? 

(4.18) A e - PAI, 

with domain 
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< 4 - 19 > P(A r ) = {(n 1 ,n 2 ,* 2 )ez R | 4 > 2 ew 1>2 , n 2 = <f> 2 (0)}. 

In [10], Burns and Cliff stated but did not prove the following theorem. We 
include the proof for completeness. 

THEOREM 4*2: The operator A R is the generator of a C^ -semigroup 

{ T R (t) lt>0 — Z R , and T R (t) = P T(t)I. 

Proof : Let I R be the identity operator on Z R . It is easy to show 

that P and I are bounded and it follows that PT(t)I is strongly 

continuous in t. Also, since p is a left inverse for I , we have that 

PT(0) I = I R . In order to show that P T(t)I satisfies the semigroup property, 
first note that the solution to (4.13) (with u = 0) is independent of the 
initial choice 4> j • Hence, for any , <j>j€L 2 (-r,0; IP), the equality 

(4.20) T(t)(n 1 ,n 2 ,<i> 1 ,<f) 2 ) = T(t)(n 1 ,n 2 ,4> 1 »<[> 2 )» t >. 0 

is obtained. Let (w(t),y(t),w t ,y t ) denote the solution to (4.13) 

a 

corresponding to (n ^ >n 2 >n^ * 4> 2 ) and u = 0 (recall with our notation x(t) = 
col(w(t),y(t))). For tj and t 2 > 0, 

.pT(t 1 )ipT(t 2 )Kn 1 ,n 2 ,+ 2 ) = PT(t 1 )iPT(t 2 )(n 1 ,n 2 ,i 1 ,+ 2 ) 

= PT(t 1 )l (w(t 2 ),y(t 2 ),y t ^) 

— PT( t ^ ) (w( t 2 ) , y( 1 2 ) , w( 1 2 ) ,y^_ ). 


(4.21) 



-33- 


In view of (4.20), (4.21) becomes 


PT(t 1 )IPT(t 2 )I(n 1 ,n 2 ,* 2 ) = PT(t 1 )(w(t 2 ),y(t 2 ),w t ,y t ) 


P(w(t + t ),y(t + t ),w ,y ) 

14 t 1 Tt 2 12 


(“ (t l + + 


A straightforward calculation yields 


PT(t 1 + t 2 )T (nj *n 2 » 4>2 ^ 


(" (t i + t 2 ).y(h + t 2 ),y ti+t2 ), 


and the semigroup property is satisfied. Using standard agruments one can 
show that the domain of the infinitesimal generator of PT(t)I is given by 
(4.19). For each (n^ .(^(O) ,<j> 2 ) e P(A R ) •*- t f°H° ws that 

PT(t)I - I 

lim (rij ,4*2(0) , <f> 2 ) = (Aq(w(0) ,y(0)) + A 1 (w(-r),y(-r)),y t ). 

t-^ 00 

It is easily verified that 


A^(nj > 4>2 ( 0 ) > ( t ) 2^ = (A Q (w(0),y(0)) + A 1 (w(-r),y(-r)),y t ) 

for (n 1 ,<t>2^ 0 ^ therefore A R is the infinitesimal generator of 

PT(t)I and T R (t) = PT(t)I. 
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The ’’reduced" abstract evolution equation becomes 


(4.22) 


z R (t) = T R (t)(n 1 ,n 2 ,* 2 ) + / T R (t - s)^ u(s)ds 


where 8 r u - (Bq u,0) (here Bg uOf 1 should be written as an element of 

1? x 1$). Theorem 4.2 and the action of the operators P and I yield the 

equivalence of (4.6) to (4.22), and hence z R (t) = (w( t) ,y(t) ,y t ) with x(t) 

= col(w(t),y(t)). The cost functional for the control problem becomes (4.7) 

with the zero operators in Q adjusted to the appropriate spaces. As a 

result of these observations, we conclude that solving the optimal control 

problem associated with (4.20) is equivalent to solving the optimal control 

problem of the original system (4.14). 

Define the reduced AVE operators by = PA N I. Computations 

R R 

similar to those in the proof of Theorem 4.2 reveal that the semigroup T N (t) 

R 

generated by Ajj satisfies T^J(t) = PT^(t)I and furthermore, it is easily 
verified that A ^ and T^(t) are the same operators that result if the AVE 
scheme is applied directly to (4.21). For (r^ , n 2 , ^ 2 >eZ R , it follows that 


T^(t)(n 1 ,n 2 > ( f > 2) = ^ tN » n 2 > <t> 2 ^ 


= P(T N (t)(n 1 ,n 2 ,n 1 »<j) 2 )). 


As commented earlier, Banks and Burns [3] proved that T N (t)z T(t)z 
uniformly in t on compact intervals for each z€Z; thus T N (t)(rij ,n 2 .tij , <t> 2 ) 

* T (t : )(n^ • 1 ^2* r ^l ^ nce P is bounded, T R (t) ■> T R (t) strongly, and 

uniformly in t on compact intervals. These results are summarized in the 

following theorem. 
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THEOREM 4.3c Let z^J(t) denote the solution of the reduced AVE scheme 

approximation to (4.22) and suppose that z R (t) = (w(t) ,y(t) ,y t ) where x(t) 

= col(w(t),y(t)) is the solution to (4.1) with initial data (n^ <f> 2 ) • 

N 

Then z R (t) * z R (t) iB. Z R , and the convergence is uniform in t on compact 
intervals. 


If the AVE scheme is applied to the reduced equation (4.20) and the 
appropriate basis chosen, then the system of equations 

[w R (t>]=A£ w£( t ) + B£ u(t) 

W R (°) = col(n 1 ,n 2 ,( ( {) 2 )^,...,( ( }, 2 )^) 

is obtained. The (n + qN) x (n + qN) matrix has the form 
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where I is th.e q x q identity matrix, 0 is the q x p zero matrix, and 
0 represents the zero matrix of appropriate size. The reduced scheme has 
effectively discarded the zero columns of A^ that appear in (4.16). 

The strength of the Chandrasekhar algorithm lies in a comparison between 
the number of differential equations in its approximation to the number of 
differential equations in the approximation of the Riccati operator. We now 
demonstrate that the Chandrasekhar algorithm drastically reduces the number of 
equations that it is necessary to solve in order to approximate the gain 
operator. Moreover, if the structure of the hereditary system is such that 
the additional reduction technique discussed above applies, then further 
computational reductions occur. 

Let the rank of the matrix Q in (4.2) be p Q . For the AVE scheme, the 
approximations Q N also have rank p Q and it is possible to obtain a 
factorization Q N = [V**] 1 V N where V N is a p Q x n(N + 1) matrix. Thus 
the set of Chandrasekhar equations associated with (4.15) and (4.17) contains 

C(N) = (m + p 0 )[n(N + 1)] 


equations. The Riccati differential equation for (4.15), (4.17) has the form 
(see Gibson [25]) 


[At)] = -Q N - [A N ] T P N (t) - P N (t)A N 
+ P N (t)B N [R] -1 [B N ] T P N (t) 

P N (t f ) = 0. 


(4.23) 
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Taking into account the symmetry of the the system, (4.23) contains 

R(N) = 

differential equations. Similar counts for the reduced system yields 

C r (N) = (m + P Q )(n + qN) 

and 

R r (N) = l l H>(n 2 1 ' + 1} , 

respectively. 

The special form of Q N yields p Q < n and thus p Q is independent of 
N. Tables (4.1) and (4.2) give a comparison of the sizes of the different 
systems for various values of the parameters. The parameters in Table (4.2) 
are for the wind tunnel problem of Example 5.2. 



Table 4.1 ? n ,— 2, q — 1, m~l, p^ — 2. 


N 

R(N) 

C(N) 

r r (n) 

C r (N) 

16 

595 

102 

171 

54 

32 

2775 

198 

595 

102 

64 

8515 

390 

2211 

198 
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Table 4.2: n=3, q=l, ra=l, p Q =l, 


N 

R(N) 

C(N) 

R r (N) 

C r (N) 

16 

1326 

102 

190 

38 

32 

4550 

198 

630 

70 

64 

19,110 

390 

2278 

134 


Based on an equations count only, for sufficiently large N (which is 
small for most problems) the Chandrasekhar equations should offer substantial 
savings in time necessary to compute the gain. A comparison of R(N) to 
C(N) shows a reduction in the number of equations from O(N^) to Q(N). If 
in addition the F-reduction technique applies, then the number of equations 
is further reduced by a factor approximately equal to q/n. 

A few remarks concerning generalizations of the problem discussed above 
are worthwhile. The inclusion of multiple discrete delays in equation (4.13) 
does not reduce the effectiveness of the Chandrasekhar equations. These 
additions affect the size and form of A N given by (4.16) (see [3]). In 
particular, for 2 delays A N becomes an n[2N + l] x n[2N + 1] matrix, and 
the equation counts above are correct if C(N) is replaced by C(2N), R(N) 
by R(2N), etc. The presence of a continuous delay term does not affect these 
equation counts; it only affects the form of A N . A numerical example 
involving two delays is given in Section V. 

Another generalization would be to include a term of the form <Gx(tj), 
x(tf)> in the cost functional (4.2), which penalizes the final state. The 
addition of this term is motivated by the fact that in the finite dimensional 
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case the inclusion of such a terra has virtually no effect on the derivation of 
the Chandrasekhar equations other than to change the final values of K(t) 
and L(t) (see [14], [31]). In this case, the Chandrasekhar equations for 

the AVE approximation scheme may be written 

K N (t) = R _1 [B N ] T [L N (t)] T CL N (t) 

K N (t f ) =R- 1 [B N ] T G N 


L N (t) =-L N (t)[A N -B N K N (t)] 


L N (t f ) = 


where 


g n = 


G 0 

0 0 


0 being the appropriate zero matrix and B^ and C are matrices satisfying 


I" tjN i T --N _N r .N iT N N a N ^ N „N -l rT ,N,T _,N 
[BjJ CB 1 = -Q - [A J G - G A+G B R [B] G. 


We point out that in the case where only discrete delays appear, the special 

structures of A 1 ^, B^, G N , and Q N allow us to find C independent of N. In 

N 

this case the matrices B^ have a constant rank Pq with p^ < (k + l)n 

where k is the number of discrete delays. In the case of a continuous 

delay, it follows that for each N, p Q = p 0 (N) and C = C N are functions 
of N. However, it still follows that Pq(N) _< (k + l)n for all N. 
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No claim is made concerning the convergence of these equations to an 
infinite dimensional form, but we do point out that the equations appear to 
converge numerically (see Examples 5.1 and 5.5), and it is interesting to note 
that an infinite dimensional version for this case has been derived in [27]. 

The use of the Chandrasekhar equations is not restricted to the AVE 
scheme. They may be applied to other schemes as well. Approximations using 
splines have become increasingly popular ([4], [6], and [7]). The spline 
scheme introduced by Banks and Kappel [6] does not satisfy the the strong 
convergence criteria on the adjoint semigroups that our work requires. 
However, a modification of this scheme by Kappel and Salamon [30] does. In 
fact, Kappel and Salamon have now applied a Chandrasekhar algorithm to their 
numerical examples and obtained satisfactory results. Another scheme for 
which the Chandrasekhar equations are applicable is the Legendre-tau 
approximation scheme introduced by Ito and Teglas in [28] and later applied to 
hereditary control systems in [29]. This scheme satisfies the strong 
convergence criteria on the adjoint semigroups and retains the low rank 
condition on Q N which makes the implementation of the Chandrasekhar 
equations effective. 


V. NUMERICAL EXAMPLES 

In this section we present several examples which illustrate the numerical 
efficiency of the Chandrasekhar equations and the F-reduction technique when 
applied to delay-differential systems. These examples are presented to show 
the efficacy of the methods and are not intended to be a complete numerical 
test. The "applications" in this section are twofold. First, the two 
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computational reduction techniques are applied to the optimal control problem 
on a finite interval. Secondly, since the infinite interval problem is the 
"limiting solution" of the finite interval problem under appropriate 
conditions, the Chandrasekhar approximations are integrated backward in order 
.to obtain a "steady-state" value. This steady-state value is then used in the 
forward integration of the state equations, and the optimal control is thus 
approximated. 

Each of the examples presented here was numerically solved on several 
different computers. In order to make a consistent comparison between 
examples, all of the computer run times recorded in this section will refer to 
the CPU time required on an IBM 3081 (located at VPI & SU). However, the 
plotted data given here may be from different machines, specifically, either a 
VAX 11/750 or the Cyber VPS-32 (Cyber 205 with enhanced memory located at NASA 
Langley Research Center) running in scalar mode. A standard fourth-order 
Runge-Kutta scheme was implemented to solve all differential equations, and 
unless stated otherwise, a fixed stepsize of h = .01 was chosen. For the 
finite interval optimal control problems, the Chandrasekhar equations were 
integrated backward from t^ and the gain values stored in increments of 
10*h. In the forward integration of the state equations, linear interpolation 
was used to obtain the intermediate gain values. 

Example 5.1 : This problem has been numerically solved in [3] and [5] and 
an analytical solution given in [3], The optimal control problem is to 
minimize 

2 2 2 

J(u) = 10x,(2) + / u (s)ds 

0 


where 
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also included an additional computation to obtain which the present 

results do not. The total CPU time for our computations was only 4.18 
seconds, and even after taking into account the differences between the two 
computers, the reductions are still quite substantial. The reduction in CPU 
time may be credited to the reduction in the number of differential equations 
solved. In [3], the authors obtained the gain via a Riccati equation and thus 
had to solve 4851 differential equations in addition to 98 state equations. 
Using the techniques discussed in this paper, it is necessary to solve only 
150 differential equations to obtain the gain and 50 additional equations to 
obtain the optimal trajectories. 

Figures 5.1 - 5.3 illustrate the convergence of the approximation 
[u N (t)] to the optimal control u*(t). Note that in Figure 5.3 the 
approximation [u^^(t)] follows the true solution closely until a time 
value of approximately 1. and then appears to lose accuracy. Implementing a 
smaller stepsize in the Runge- Kutta scheme corrected this. 

Also included here are graphs of the approximations to the functional 
gain a + iqC^a), of the gain operator K(t). Recall that the averaging 
scheme produces piecewise constant approximations to functions in L, 2 (-r, 0 ); 
thus the elements O and A in Figures 5.4 - 5.11 represent these constant 
values on the appropriate subintervals [a^_^,a^] I - r ,0 ] * In order to 

obtain illustrative clarity, these points were connected to form a smooth 
curve for N = 100. Pictured are approximations of the functional gains 
K^(t,a), for N = 8, 32, 100, at the values of time t = 0, .25, .5, .75, 1., 
1.25, 1.5, and 1.75. These figures show not only the convergence in N of 
the functional gains, but also the evolution in time of these approximations. 
Note that the approximate functional gains calculated at t = 1.25, 1.5, and 
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1.75 do not clearly Indicate the discontinuities that actually appear in the 
functional gain due to the presence of a terminal cost terra in the functional 
to be minimized (see Delfour [20]). This phenomenon will be discussed in more 
detail in Example 5.5. 

Finally, graphs of the resulting state approximations for N = 8 and N = 
100 which show rapid convergence of the states are given in Figures 5.12 and 
5.13. 


Example 5.2 : This next example is an application of the Chandrasekhar 
equations and of the F-reduction technique to an infinite time optimal 
control problem. This illustrates that the Chandrasekhar equations may be a 
viable method of computing the constant gain associated with the infinite time 
problem. The motivation for studying this particular example is not only that 
the problem has the special structure which fits our framework, but that it is 
a practical problem as well. The system is a model for fine tuning the raach 
number in a cryogenic windtunnel (National Transonic Facility) constructed by 
NASA at its Langley Research Center in Hampton, Virginia. This model has been 
studied in [37], [39], and the LQR problem numerically approximated in [7], 
[33]. 

The controller of the system is an actuator attached to a wind guide vane 
and finely tunes the raach number by changing the angle of the vane. The state 
consists of the variation of the raach number, the variation in the guide vane 
angle velocity, and the variation of the guide vane angle. The equation has 


the form 
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-a 

x(t) = 0 

0 

x(s) = col(-. 1,0. ,8.547), -.33 < s < 0, 

2 

where the parameters 1/a, m , k, and b have the values 1.964, 36., -.0117, 
and .8, respectively. The quantity to be minimized is 

00 

J(u) = J (x^(s)Qx(s) + u^(s))ds, 

0 

where Q = diag(10^,0,0). 

In [7], the authors were comparing the AVE approximation scheme to a 

spline approximation scheme. They applied a Newton iteration scheme to the 

Riccati matrix equations that resulted when the AVE scheme was used (for N = 

2, 4, and 8). However, the Newton iteration scheme did not converge to the 

solution of the Riccati equation for N = 8. 

In our approach, we solved the finite time optimal control problem on an 

interval [0,tf] of sufficient length so that the gain K N (t) satisfied 
N N + 

K (t) ->■ Kq as t + 0 , (i.e., we integrated backward to a steady-state 
N 

solution Kq). This method resulted in a convergent scheme for N = 8, and in 
fact, for larger values of N as well. 

48 

Figure 5.14 illustrates the convergence of the functional gain (0,a) 
to a steady-state value for the values of tf equal to 2, 3, 5, and 10. Note 
the rapid convergence of the gain to a steady-state value and that the values 
for tf = 5 and 10 are almost indistinguishable. 
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The states and optimal control computed by using the steady-state gain, 
48 

K (0) are given in Figures 5.15 — 5.18. As can be seen, the response of the 
system is very good, and the states are driven to zero rapidly. These results 
agree well (pictorially) with those found in [7]. 

The total computer time required to compute the gain, optimal control, and 
states for this problem (N = 48) was 7.76 seconds. This computation used 
both the Chandrasekhar equations and the additional reduction technique 
outlined earlier. The same problem was also solved on the computer using only 
the Chandrasekhar equations and not the additional reduction. In this case 
the total CPU time was 21.55 seconds, a significant increase. 


Example 5.3 : In this example the Chandrasekhar algorithm and F-reduction 
techniques are applied to the hereditary model of the two-dimensional airfoil 
discussed in Section 2. The resulting delay model has the form (see equations 
(2.16) and (2.17)) 


(5.1) 

with output 

(5.2) 


x(t) = A q x(t) + A 1 x(t - s) + Bu(t) 
x(s) = <j>(s) -r _< s _< 0 


y(t) = Cx(t). 


Here, equation (2.17) has been multiplied by the inverse of E. The matrices 
in (5.1) - (5.2) are 
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“-7.3865 

-33.8517 

-4479.50 

-5711.27 

3.0804 


.9378 

-12.4685 

568.73 

-7068.77 

-.3911 

> 

o 

ii 

1.0 

0.0 

0.0 

0.0 

0.0 


0.0 

1.0 

0.0 

0.0 

0.0 


-2.7443 

1404.29 

-1664.26 

-40701.67 

-253.856 


Aj (5,5) = -47.00, A^i.j) = 0, (i,j) * (5,5) 

B = [-82.333 191.286 0.0 0.0 864.533 ] T 

C = diag[/50.00 /5070 /ToTO /TO l] 

<|>(s) = col[-.80 .50 .055 .029 50.0], se[-r,0], 

and the time delay is r = *05. The cost functional that we desire to 

minimize is 

t f 

J(u) = / (y T (t)y(t) + Ru 2 (t))dt 

0 

where the final time is t^ = .25 sec., and R = 10. 

Because of the small delay, the stepsize h used in the Runge-Kutta 
scheme was reduced to .001. The computations for this particular example 
exhibited in Figures 5.19 - 5.24 were performed on a VAX 11/750 computer using 
single precision arithmetic. 

The parameters used to construct the above matrices are a slight variation 
of those found in [10]. In particular, the non-diraensionalized distance, x^, 
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from the airfoil center of mass to the elastic axis (e.a.) is .5 which places 
the e.a. at the quarter— chord point. In addition, the free— stream velocity 
U was chosen to be 1500 in/sec. This particular value was selected since it 
produced unbounded oscillations in our open loop system (a few numerical tests 
determined a "flutter speed" of approximately 1375 in/sec for this model). 
This unstable plant models the airfoil operating under flutter conditions, 
i.e., the pitch and plunge modes couple in a manner that allows unstable 
oscillations to occur. Since it is desirable to eliminate flutter, the 
controller should be able to damp effectively the pitching and plunging 
motions. 

The reduced AVE scheme with N = 16 was implemented in order to 

approximate (5.1) - (5.2). The time-varying gain values produced by the 
Chandrasekhar equations were used to obtain the optimal control, u 16 (t). 
Graphs plotting the closed loop responses (dotted graph) of the LQR design are 
compared with the responses of the open loop (solid graph) system in Figures 
— 5.23. As can be seen in the figures, the pitching and plunging 
motions (x 3 and x^ respectively) and their velocities (xj and x 2 ) are 
quickly driven to zero in the closed loop response, illustrating the 
effectiveness of an LQR approach for our model. The optimal control which 
produced these responses is given in Figure 5.24. 

The computational advantages of the reduction techniques discussed in this 
paper are very apparent in this example. For N = 16, in order to solve for 
the gain via the Riccati equation, one must solve 3,655 differential equations 
(this includes taking advantage of the symmetry of the equation) ; moreover, 
there are 85 state equations to solve. In contrast, if the two techniques 
discussed here are employed, it is necessary to solve only 126 equations to 
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obtaia the gain and 21 equations to obtain the state. If N is further 
increased to 48, then the number of equations becomes 30,135 Riccati and 245 
state equations compared to 318 Chandrasekhar and 53 state equations. A 
similar problem (i.e., same system with different aerodynamic parameters) was 
solved on an IBM 3081. For N = 16, the total CPU time required by the 
methods here was 8.47 seconds. When N was increased to 48, the CPU time 
increased to 20.13 seconds. 

Example 5.4 ; This example is a variation on Example 5.3. The parameters 
x^ and U defined previously are set equal to .297 and 1325., 

respectively. This value of x^ moves the e.a. from the quarter chord point 
(note also that from the figures U is below the "flutter speed"). The 
resulting matrices for equations (5.1) - (5.2) may be found in [12]; all other 
values are as in Example 5.3. The purpose of this example is to compare the 
state response obtained by using the steady-state value of the gain to the 
state response obtained by using the time-varying gain. In applications, use 
of a constant gain is often preferable since the gain values do not have to be 
continually updated at each instant of time. 

The Chandrasekhar equations were integrated on the interval [0, .25] , and 
the value of K^(t) at t = 0 was used as the steady-state value. Figures 
5.25 - 5.29 compare the open loop response to the closed loop response on 
[ 0 , . 5 ] obtained by integrating the state equations using the constant gain. 
The closed loop responses produced by the time-varying gain yielded results 
that appeared graphically the same as the responses produced by the constant 
gain. There were, however, slight numerical differences. Graphs of the 
absolute values of these differences for h(t) and a(t) are given in 
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Figures 5,30 and 5.31. These graphs are representative of the results for all 

of the states. Figure 5.32 is included to show how the kernel portion of the 

48 

functional gain (t,a) evolves in time. In [1], the authors treated a 

similar problem. However, their results showed high frequency oscillations 

near t = .25, which may be due to numerical instabilties. 


Example 5.5 : This last example is a one-dimensional equation with two 

delays considered in [20]. The optimal control problem is to minimize 


J = 10x 2 (2) + J u 2 (s)ds 
0 


where x(t) satisfies 


x(t) = x(t - .5) + x(t - 1.) + u(t) 
x(s) =1, -1 < s £ 0. 

This problem illustrates the type of discontinuities that may occur in the 
kernel of the gain, and how they are affected by the approximation scheme. 

Since the problem has delays of 1. and .5, Kj(t,a) will have a fixed 
discontinuity at a = —.5 for each t. However, because of the presence of a 
final cost term, Kj^tja) will also have discontinuities at (see Delfour 
[20]) aj = ~1 • + tj - t and a 2 = -.5 + t^ - t (for values of t which 

place a^ or ^ [~1 ,0] ). The approximations K^^(t,a) at t = 1.9, 

1.4, 1., and 0. are given in Figures 5.33 - 5.36. Note that Figure 5.33 

should show discontinuities at a = -.4 and a = -.9, and Figure 5.34 should 

show a discontinuity at a = —.4. The smoothing that occurs is due to the 
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averaging of values which the Runge-Kutta scheme uses to produce each new 
iterate, and is not a characteristic of the Chandrasekhar algorithm. Though 
not exhibited here, use of a lower order Euler scheme (which does not average 
values for the next iterate) shows definite jumps at the points of 
discontinuity, but sacrifices overall accuracy. The effect of the smoothing 
of the kernel is minimized, however, since the action of the kernel is through 
an integral term. Consequently, accurate values for the optimal trajectories 
and optimal control are still obtained. 
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APPENDIX 


PROOF OF LEMMA 3.1 

Let R^(t) and [8^1 (t) be the Nth approximation analogues of (3.5) 

o S 

and (3.7). Note that R N (t) is self-adjoint and satisfies llR N (t)ll > M for 

s s — 

each N, s, and t. It follows that II [R^] -1 (t) II _< £ uniformly in N, s, 

and t. Also, it follows that (T g »)(t), (Tg*)(t), and (F g «) are uniformly 

bounded in N, s, and t. Combining these observations with the uniform 

bounds on B N and Q N , it follows that llB N (t)ll is uniformly bounded in N, 

s 

A N — 1 * 

s, and t. Therefore, II [R g ] ( t ) II is also uniformly bounded. Recalling 

that [R g ] (t) is self-adjoint and that the norm of a bounded operator 
equals the norm of its adjoint, we have established the lemma. 


PROOF OF LEMMA 3.2 

We first observe that strong, pointwise convergence of T N (t) and 

[T^(t)] along with an application of the dominated convergence theorem 

implies that (T^»)(t), (T^ »)(t), and (F N *) converge strongly and pointwise 
s s v s J * 

t0 (TgOU). (T s *)(t), and (F g »), respectively. In addition, (T g «)(t), 

(T^ »)(t), and (p g »)(t) converge strongly and pointwise to (T «)(t), 

(T s .)(t), and (F g *)(t:), respectively. Since B^, and converge strongly, 

it follows that R^(t) -► R (t) and s!?*(t) •> 8*(t) strongly and pointwise 

s s s s 

for 0_<s_<t_<tjr. Taking the adjoint of (3.6) and using similar reasoning, 

*N * 

it can be shown that B (t) + B (t) strongly and pointwise for 

s s i ^ 

0 £ s £ i t f • The identity R g (t) - R~ X (t) = R g (R. - R g )R g X (t) now 

*N"”^ ''—i 

implies that R (t) + R (t) strongly and pointwise. Therefore, it now 
s s 
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*N A N 1 

follows that 8 R (t) * 8R (t) strongly and pointwise for 

s s 

0 s t t^. As indicated previously, the convergence of 
S(t,s) was shown by Gibson [25], Using equation (3.8) 
S^*(t,s) may be expressed as 

<A.l) [sVs)]* * = [At - s)f z - /' S s [R«]- 1 <n)[^]*[i N ( C 

s 

Combining the assumptions on T N *(t) and B N * with Lemma 3.1, 

of (A.l) is seen to be uniformly bounded in N, s, and t. As a 

comments above, the integrand also converges pointwise to 
* *-l, * * 

B R (n)B T (z - n)z, and an application of the dominated 

o S 

theorem yields S (t,s) + S (t,s) strongly and pointwise. 
boundedness of S^*(t,s) (and hence S^(t,s)) follows directly 

PROOF OF THEOREM 3.3 

An integration of (3.14) results in the equation 

c f 

(A.2) K N (t)z = J R~ 1 [B N ]*[L N (n)]* L N (n)zdn. 

t 

Equation (3,15) may be rewritten as 


L N (t)z = -L N (t)A N z + L N (t)B N K N (t)z 
and a variation of parameters formula yields 


S^(t,s) to 
we see that 

i* 

- n) J zdq. 

the integrand 
result of the 

convergence 
The uniform 
from (A.l). 
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(A.3) 


L N (t)z = V N T N (t, 


- t)z - / L(n)B N K N (n)T N (n - t)zd n . 


Using the assumption that L N (t) = V N S N (t f ,t), (A. 2) and (A. 3) become 


(A. 4) 
and 


K N (t)z = / R -1 B N * s N *(t n)v N * v N s N (t ,n)zdp. 


(A. 5) V N S N (t f ,t)z = V N T N (t f - t)z - / V N S N (t f ,ri)B N K N ( n )T N (n - t)zd n , 


respectively. The dominated convergence theorem, Lemma 3.2, and Gibson's 
convergence results (see [25]) for K^(t) imply that the limit as N °° of 
(A. 4) exists and 


(A. 6) 


K(t)z = 


r 1 B* s*(t f ,n)v* vs(t f ,p)zdn. 


Since K^(t) and T^(t) are each strongly continuous in t and converge 
to K(t) and T(t) in the strong sense, they are uniformly bounded in N 
and t on compact t-intervals. Therefore, another application of the 
dominated convergence theorem and Lemma 3.2 imply that the limit as N + °° Q f 
(A. 5) exists and 


(A.7) vs(t f ,t)z = VT(t f - t)z - / vs(t f ,nBK(n)T(n - t)zd n . 


Defining L(t) = VS(tf,t) we have shown that the pair K(t), L(t) satisfy 

the Chandrasekhar equations (3.16) - (3.17) and K^(t)z + K(t)z for all 

z£ff. If U is finite dimensional, then it follows (see Gibson [25]) that 
N 

K (t) * K(t) in the uniform operator topology. 
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PROOF OF THEOREM 3.4: 

A A 

Let (K(t),L(t)) be another solution pair for (3.16) - (3.17). It 
follows from (3.16) that 


II (K(t) - K(t))zll < / hr " 1 B*(L*(n)L(n) - L*(n)L(n))zHdn 

t 

"I ^ A ^ 

_< HR B U / il L (n)L(n) - L (n )l^( n ) ll iiziidn, 

t 


which implies that 


iiK(t) - k( t ) ii _< hr " 1 b »/ HL*(n)L(n) - L*(n)L(n)iid n . 


Adding and subtracting the appropriate term yields the inequality 


nK(t) - K(t) ii < hr " 1 b*hJ (HL*(n)L(n) - L*(n)L(n)n + iiL*( n )L(n) 


- L (n)L(n) ll)dn 


— 1 ^ ^ a A^ ^ ^ A 

< IIR 8 II J (iiL(n)n HL (n) - L (n) ll + llL (n)n iiL(n) 


- L(n)u)dn. 


Since the norm of an operator equals the norm of its adjoint, it follows that 


-1 * 


llK(t) - K( t ) II < HR 8 ll/ (llL(n)ll + llL(n)ll) ||L( n ) - L(n)lldn. 
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The principle of uniform boundedness implies the existence of a constant 
0 such that 


(A. 8) 


A 


llK(t) - K(t) II _< Mj / 


A 

nL(n) - L(n ) ii dn - 


Similarly, equation (3.17) yields the inequality 


ii£(t) - L(t) ii < / iL(n)BK(n) - L( n )8K(n)ii HT( n - t) iidn 

t 

tf 

< / HT(q - t)ll HL(q)BK(q) - L(q)BK(q) + L(q)BK(q) 
t 

" L(ri)BK(n)lldn 


< / HT(n - t ) II (llL(q)BII llK(n) - K( n )|| 
t 


A 

+ IIL(q) - L ( n ) II IIBK(ri) ll)dq. 


Again, the principle of uniform boundedness can be applied to obtain the 
existence of M 2 > 0 such that 

(A. 9 ) »L(t) - L(t) ii < 1^2 J (iiL(n) - L(n) ii + nK(n) - K(n) II )dn. 

t 

Substituting (A. 8) into (A. 9) yields 

t £ 

f ~ t A 

(A. 10) IIL(t) - L(t) II < M 2 f ( llL(n) - L(n) II + M { f ||L(s) - L(s)||ds)dq. 

t q 
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Since t _< n in (A. 10), the lower limit may be extended to obtain 


HL(t) - L(t) II < M 9 J ||L(n) 

t 


- L(n)ndn + M x h 2 / 


J llL(s) - L(s) lldsdn. 

t 


Applying Fubini's theorem and interchanging the order of integration, we have 


A “ A 

llL(t) - L(t)ll _< (M^ + M^t^ " / HL(s) - L(s)llds, 

t 


and hence 

(A. 11) HL(t) - L(t ) II < (M 2 + M x M 2 t f ) J HL(s) - L(s) lids. 


In order to justify the existence of each of the integrals appearing above, we 
note that each operator is strongly continuous on [0,tf], and hence the norms 
of these operators are lower se mi-continuous on 1 0 , t £ J (see Kato [34]). 
Lower semi-continuity implies measurability, and since the norms are uniformly 
bounded, they are integrable. 

Gronwall^s inequality holds for integrable functions almost everywhere 
(Reid [41]), hence (A. 11) and Gronwall's inequality imply that 
HL(t) - L( t) I! = 0 almost everywhere in [0,t f ]. Since llL(t)z - L(t)zll 
continuous for each z€H and equal to zero almost everywhere, it follows that 

A 

L(t) = L(t) everywhere on [0,tj]. Equation (A. 8) immediately yields 

A 

K(t) = K(t) for t£[0,t f ] , and the proof is complete. 
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